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ABSTRACT 



The problem of maximum range atmospheric reentry for an orbiting 
lifting glider was treated by Bryson and Denham'*’ by the method of 
"steepest ascent". The same problem is undertaken here by a method 
of differential corrections developed by Faulkner, This method makes 
use of a Newton-Raphson type iteration based on paths which satisfy 
the Euler-Lagrange equations, A comparison of results is made,, show- 
ing large differences in control variable history , and longer range 
for the path obtained by differential corrections. The problem was 
characterized by a sharp "ridge" in the domain of the starting values 
of the adjoint variables and the effect of this on the convergence of 
both methods is discussed. Finally , the difficulty of choosing initial 
approximations for the starting values of the adjoint variables is 
discussed, and a method is presented for obtaining these from a nominal 
path as the first step in the computer routine, 
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A, E, Bryson and W, F. Denham, A steepest-ascent method for solv- 
ing optimum programming problems (Ratheon report BR-I303), Raytheon 
Company, Bedford, Mass,, 1961, 
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INTRODUCTION 



A problem of optimum control may be described very roughly as 
the determination of the decision-making or control function which 
will result in the largest possible value of the final "payoff" from 
some continuing process 0 

The process might be illustrated as shown in figure l<,ls where 
P is the "payoff" , something whose value depends on the terminal point 
or the path history 8 and Q is a set of given constraints which must be 
satisfied,, 




The methods of solution commonly used now were foreseen at least 
forty years ago „ but the volume of computation required limited their 
application to only the simplest of problems „ In the past ten years 
there has been an explosion of interest due to the availability of 
high-speed digital computers on the one hands and on the other hand 
to the urgent need for solutions in such applications as space 
trajectories „ The "state of the art" is such that no method is com- 
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pletely general, and a great diversity of approaches is employed, 
reflecting the absence of a firm framework of mathematical theory 
concerning nonlinear differential equations. 

In this thesis a comparison is made between two methods for cal- 
culating optimum angle of attack programs for a lifting unpowered 
vehicle to attain maximum range, starting from circular orbital 
velocity at a height of 300 thousand feet above the surface of the 
earth. The problem was chosen because it is typical of a class of 
problems in which there is now great interest, and because two of 
the foremost workers in the field of optimization, Arthur E, Bryson 
and W, F. Denham, had published a clear and well-documented solution 
by the method of "steepest ascent". 

The comparison was suggested by Stanley Ross, while one of the 
writers was engaged in summer field work under the sponsorship of 
the Lockheed Missile and Space Company, John V, Breakwell arid George 
Leitmann generously gave many hours of their time to general dis- 
cussions and made specific suggestions which were very helpful. 
Professor Frank D, Faulkner of the U, S„ Naval Postgraduate School 
not only laid out the general form of the problem and overcame the 
difficulties that arose, but also the method of solution used was 
the one developed by him from the basic method of differentials of 
Bliss, 
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Chapter I 

GENERAL FORMULATION OF A PROBLEM OF OPTIMUM CONTROL 

1.1 The equations of motion, the control variables and the con- 
straints . 

It will be helpful to define an optimum control problem in 
somewhat more exact terms than those given in the introductory 
remarks . 

Suppose we are given a system whose behavior is described 
by a system of differential equations which we shall call the 
equations of motion. We will suppose that wherever any of these 
differential equations is of higher than first order, we have de- 
fined such additional variables as necessary so as to obtain n 
first order differential equations of the form 

= S2* s n j, *0 

i ™ 1, 2, o.o, n 

where the variables s. are the state variables, the p^ are the 
control variables, t is the independent variable, and the dot 
denotes differentiation with respect to the independent variable. 
It is assumed that the f^ are of class C" . 

It will be convenient for what is to follow, if matrix 
notation is adopted at this point. The equations of motion may 
be written as the single matrix equation 
(1.1) S = F(S , P, t) 
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where 



S = . , an nxl matrix of the dependent state 

variables 



n 



1 



F = . 




an nxl matrix of known functions of the 
control variables, the state variables, 
and possibly the independent variable 



P 



' m 



an mxl matrix of control variables „ 
which we are free to choose as functions 
of t. 



In addition, there are k constraints or end 
e^(S,t) = 0, which may define values of the 
end point t = T, or may be of integral form 



conditions, of the form 
state variables at the 
such as 




dt = constant. 



In the latter case, we shall define an additional variable s n +]_ 



so that 




with the initial and end conditions 



4 



( 0 ) = 0 



s 



n+1 



s, 



'n+1 (T) = the given constant 



These may be written as 



e, 



1 



(1.2) E = 



, a kxl matrix of constraint Tun ct ions 



We wish to discover the particular control variable 
program which will drive this system from a given initial point 
to some terminal point which, of the set of all such points which 
satisfy the terminal constraints, gives a maximum value to some 
function, M(S,t), of the state variables; that is, for the ter- 
minal point t = T, possibly unknown, 



1.2 Bounded control variables and inequality constraints. 

The permissible domain of the control variables may likewise 
be limited. Some examples of restrictions on the control for com- 
mon physical systems might be: 



M(S,t).j> = M(T) = max. 



^ ^ P max. 
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P = P max for 0 < t < t 1 and 
P = 0 for t-j<t< T. 

Bounded control variables appear frequently in real physical 
systems. For solving these problems it may be useful to re- 
place the inequalities by mathematically equivalent equalities. 
This is accomplished by the introduction of suitable new vari- 
ables in the manner due to Valentine (ref.l). For instance, 
for a constraint 

0 £p£P max 

one may define the real variable g, so that 

P( p max - P) - g 2 = 0« 

This procedure transforms the inequality to an equality, and we 
adjoin this equation to the set (1.1). The same procedure is 
applied in the case of inequality constraints or end conditions. 

We will generally require that all state variables be speci- 
fied initially, but this is not necessary. One or more of these 
may be unspecified, that is "free". 

The general form of these problems, then, is the two-point 
boundary value problem. With the stipulation that the functional 
M is to be maximized (or minimized) they become, depending on the 
formulation, the classical problem of Bolza or problem of Mayer 
of the calculus of variations. Since to maximize some quantity, 
say u, is equivalent to minimizing the quantity minus u, we will 
henceforth consider that the term "maximum" includes both cases 
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unless otherwise specified. 

1.3 Degenerate problems. 

Some important observations may be made here. First, as 
pointed out by Faulkner (ref. 2) and Breakwell (ref. 3), the case 
in which all of the functions f^ involve any function of a par- 
ticular control variable only linearly are said to be degenerate . 
These degenerate problems typically are of the "bang-bang" type 
in which the control changes dis continuously. The neglect of 
induced drag in aerodynamic problems for instance often results 
in bang-bang control. 

Second, the constraints must not involve the control variables, 
or the problem is overspecified. If the control program is to be 
such as to produce an extreme of the functional M, then it may not 
simultaneously satisfy any prescribed constraint at any isolated 
point of the path, else the solution is generally discontinuous. 

1.4 The adjoint variables. 

We shall introduce a set of Lagrange multipliers, as yet un- 
determined 




a lxn matrix of adjoint variables. The term •adjoint" is used, 
since they will be chosen to be solutions to the system of differ- 
ential equations which is adjoint to equations for the variations 
of the given system. 

Suppose we write (1.1) as 
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S - F = 0 

and integrate its product with the adjoint vector between the 
limits of the independent variable. For convenience, we will 
call these limits zero and T, and we have, 

,T * 

/ U(S - F) dt = 0. 
o 

The corresponding variational equation is 

/ T U(«*S - F dS - F f?) dt = 0. 
o s p 

Now this is integrated by parts to eliminate from the integrand 
the variations of the state variables. 



(1.3) |u£sl T = f T f( U+UF )SS +U F 6 p] dt-fcj fI tk+ dt 

L. _J o ; 0 L S p • C K -J t k- 

where F s and Fp are the matrices of partial derivatives. 



(1.4) 



F e = 






• o 

9t n /d s x . . . af n /as n 



, an nxn matrix. 



f p = 



^l/3Pl * • * 
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• • • af i/ap m 



, an nxm matrix, 



and t is a symbol for any point or set of points where S or U 

K 

is discontinuous. 

Now, if U is chosen as a solution to the equation 
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(1.5) 



U + U F = 0, 
s 



then ( 3 ) becomes 




k 



This sequence of operations defines the system which is adjoint 



Hereafter, no distinction will be drawn between a vector in 
k-space and a 1 x k or kxl matrix, or a corresponding 
column or row of a matrix. 

1.5 Green's formula and the Euler-Lagrange equations. 

Equation (1.6) is the fundamental formula relating the varia- 
tions of the end values of the state variables with the variations 
of the control variables. Note that the formula shows the ter- 
minal value of the adjoint variable to be the sensitivity of the 
terminal value of the corresponding state variable to a change in 
control. This interpretation will be discussed further in con- 
nection with the transversal conditions (sect. 1. 9 ). 

Equation (1.6) is often called Green's formula (one form). 

See Coddington and Levinson (ref. 4) page 86. The last term on 



to the variational equations of ( 1 . 1 ), Ss-F <$S = f6p , 

s p 



and forms the 1 x n matrix, or vector 
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the right applies only at points of discontinuity of F, and 

o 

since these are discontinuities in S, they are •comers". 

Hereafter, it will be assumed that the solution does not have 

a comer unless it is otherwise mentioned. In this case the 

last term drops out. The fundamental lemma of the calculus of 

variations states that, given a control program P(t) such that 

the constraints E are satisfied, then if the coefficient of &P 

in (1.6) does not vanish for some solution U to the adjoint 

it is possible to construct a variation such that the constraints 

are still satisfied, but the function M exceeds the value on the 

first path. A proof of this lemma is included in refs. (6) and (10). 

Thus for M to have at least a stationary value M* , 

(1.7) UF = 0. (The Euler equation) 

P 

Equation (1.5) defines the system of equations which is adjoint 
to the variations of the equations of motion, (1,1). The adjoint 
equations and equation (1.7) will be called the Euler-Lagrange 
equations. It will be shown that these, together with the con- 
straints and the equations resulting from the transversal condi- 
tions comprise a complete set which determines the optimum solution. 

1.6 The maximum principle of Pontryagin. 

The functions f^ are often composed of terms which are func- 
tions of the state variables alone and terms which contain also 
the control variables. Suppose we collect terms of the first type 
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into an n x 1 matrix we will call V, and those of the second type 
into an n x 1 matrix G. Then (1.1) becomes 



S = V + G. 



Now consider these matrices to be vectors in an n-dimensional 
hyperspace of the state variables, and a particularly fruitful 
geometric interpretation of (1.7) appears. Let us consider equa- 
tion (1.7) 

U • F p = 0 

as the condition that 

U * G = max 
P 



that is, let us require that the vector G be everywhere so dir- 
ected that the projection of G upon U is maximized. The control 
program which maintains this relationship (while satisfying the 
constraints) is optimum at least locally. This statement of the 
condition is due to the Russian mathematician Pontryagin, and it 
leads to particularly enlightening geometric formulations in many 
problems. Note that the shape of the domain of G then gives in- 
sight into the behavior of the control program. 

In figure 1.2 a hypothetical domain of G is shown by the 
dashed line. The scalar quantity U • G is maximized when G has 
a maximum projection on U at all times. Figure 1.3 shows the 
situation where more than one G vector satisfies the maximum 

condition. This generally yields a discontinuity of the control, 
hence a "corner". 
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The functions g^ will be called the "forcing functions". 

It is the vector whose components are these functions which must 
be directed so as to maximize its projection on the adjoint vec- 
tor. The equation or equations which express this condition in 
the form 

Pi = h(U, S, P, t) 

are obtained try ordinary calculus from (1.7). They are sometimes 
called "steering equations". 

1.7 The Hamiltonian. 

From the foregoing, one sees that the product UF has special 
significance when U is chosen according to (1.5). For convenience, 
we shall define this quantity as the Hamiltonian function, H. 



H = UF. 



The maximum principle becomes 

(1.8) H = max , 

P 

and the general condition for an extreme of H is 

(1.9) H = 0, 

P 

which is equivalent to (1.7). 

For any one extremal, if t does not appear explicitly in the 
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equations of motion, the Hamiltonian is constant. This may be 
shown readily, for 



H = UF + UF 



= UF + UF s S + UFpP + UF t 
= (U + UF )F + UF_P + UF. . 

S P K/ 



But U + UF = 0 by equation (1.5) » and UF = 0 by (1.7), hence 

s p 



p 



(1.10) H = UF t = H t 



In a great many problems, the independent variable does not 
appear explicitly in the equations of motion, and the Hamiltonian 
becomes a constant of the system. In the next section it will be 
shown that in this case if the final value of the independent vari- 
able , T, is free (that is, it is not the quantity to be optimized 
and is not constrained) then H = 0 everywhere on a given path, 
provided the Euler-Lagrange equations (1.5) and (1.8) are every- 
where satisfied on the path. 

1.8 Systems which are linear in the state variables. 

Consider the functional 



If the differential equations of motion are linear in the state 
variables, then the functional does not involve the state variables. 
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and an absolute maximum may be obtained. If the differential 
equations are not linear, then there are no corresponding gen- 
eral proofs, ana we can look only for a stationary value and a 
local extreme for the functional. 

A system of differential equations which is linear in the 
state variables with constant coefficients may be written as 

S = A S + Q 

where S is the nxl matrix of equation (l.l), Q is an nxl matrix 
whose ith component is the forcing function of the control vari- 
able in the ith direction, and A is an nxn matrix of constants. 
Now, in analogy with the derivation of equation (1.6), we may form 
the integral 

f U(S - AS - Q)dt. 

O 

Integrating by parts 

_ T T 

jus] = / Eu + UA )S + UQ"| dt. 

0 0 

o 

Taking U to be the solution to the adjoint equation U +UA = 0, we 
have the form of the Green-Lagrange formula 

T T 

(1.11) [USJ = /UQdt. 

o o 

Suppose that we have by some means discovered a set of 
initial values for the adjoint variables and integrated forward 
until, at some time T, 
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H 1. the curve satisfies the constraints on 

the state variables and the control vari- 
able, 

H 2. the forcing function Q is such that UQ is 
maximum everywhere on the curve, and 
H 3. the terminal value of the kth component 
of the adjoint vector is unity, and all 
other components are zero at t = T, 

In this case, Faulkner (ref. 5) shows that the curve furnishes an 
absolute maximum for the kth state variable under the given initial 
values and the constraints. Akheizer (ref, 6, sect. 2.15) and 
Edelbaum (ref. 7» sects. 1.1 and 1.2) discuss conditions for 
strong and weak extremes in this and in nonlinear cases. 



1.9 The transversal condition. 

The extremals have been discussed, as the solution curves 
which satisfy the system 

C 



( 1 . 12 ) 



S - F = 0 

< U + UF = 0 
s 






H 



P 




= o 



Consider the system already discussed, and the corresponding 
Green-Lagrange variational equation (1.6), 



_ T T 
U6sl = / (UF 6P)dt 
- “* o o P 




t k + r t 
dt, =J OH dt 



V 



k & 
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We require that the tern on the left, evaluated at the starting 
point be zero. If s^(o) is specified, then 6s^(0) =0, If s^ 
is not specified, then we shall choose the particular solution 
to the adjoint so that U^(0) = 0. It is necessary that the last 
tern on the right be zero, else changing dt will offer a better 
trajectory. There are variations for which c>P is not "small" and 

r T 

J 8h dt <0, 

O 

but for first order effects, the only condition on £>S for an ex- 
tremal is that 

(U 6 S) = 0. 

The form of the equation does not depend on the constraints „ hence 
we will define as admissible paths those which satisfy the constraints. 
Since the constraints must be satisfied at t = T, the endpoint S(T) 
must lie on the n - k dimensional surface 

E(S,T) = 0 

Hence, for an admissible differential at the endpoint, 



de. = 
x 






Be 



i /as j 



ds 



3e./d tdT = 0 

i J T 



or in matrix notation, 



1 ? 



( 1 . 13 ) 



dE = (EgdS + E t dT) T = 0 

where dE and are of dimension kxl, S s is kxn and dS is nxl. 
E s may be considered to be the gradient of E in n-space. If 
E s is bordered on the right by the column E t , this may similar- 
ly be considered to be the gradient of E in the n+1 space (of 
the state variables and t) which we will call VE* . 



VE* = 



ae 1 /as 1 . . .ae 1 /as n ae^/at 



• 0 9 

de^/as^ • • • de^/ ds^ ae^/at 



Now if dS* is taken to be 



dS* = 



ds. 



ds. 

dT 



n 



( 1 . 13 ) may be rewritten as 



dE = (VE*dS*) T = 0. 



In the same manner, we may have the gradient of M in n+1 space, 



VM* = 



2K/asj 3M/ds 2 , . . _aii/ds. 



n dM/atJ 



and the augmented adjoint vector 

U* 



" \^1 U 2 * * * \ ( " H ^] 
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The vector U* must lie in the manifold spanned by the k+1 vectors 
Ve^* VM , hence it may be expressed 

= C Z* 

where C is a lx (k+1) row vector of constants and Z* is the 
(k+l) x (n+l) matrix obtained by bordering VE* below with V M* . 
Since M is assumed to be independent, Z* has rank k+1. The 
transversal condition may be stated as the condition that the 
matrices V E* , Z* and I* 

, a (k+2) x (n+l) matrix 

all have rank k+1 . 

As a simple example, suppose we have the state variables 
w, x, y and z, with all initial conditions specified and for 
t=*T w constrained, x to be maximized, y, T and z free. We may 
write 

w x y z T 
10000 constraints 
01000 maximum 
U 1 u 2 u 3 ^ adjoint variables 




U* = C 



VE* 

VM* 
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a 3x5 matrix which must be of rank 2 at t = T. Hence 

u 3 = \ = -H = 0 at t = T. 

In the development above, free use was made of material from 

11 

a research paper of Faulkner, '.'/hen each e^ and M involves only 
one variable of the set (S,t), the condition may be summarized as 
follows : 



(1.14) 



s^(T) free implies u^(T) = 0 
T free implies H(T) = 0 
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Chapter II 

SOLUTION OF A GENERAL OPTIMUM PROGRAMMING PROBLEM USING 
DIFFERENTIAL CORRECTIONS 



2.1 The missing constants of the adjoint set. 

It has been shown that if relations (1.12) and (1.14) are satis- 
fied, then the resulting path generally furnishes an optimum to the 
state S(T). Finding this particular path depends upon finding a 
particular set of constants, which are the initial values of the 
unknown adjoint variables. This being true, suppose for the moment 
that the means exist for correcting a set which is "pretty good 1 ', 
until the final "perfect" set is obtained. That is, we visualize 
in the domain of U(0),the adjoint vector evaluated at t=0, the 
point U*(0) whose coordinates are these "perfect" starting values. 

If U(0) = U*, then integration of the first two of equations (1.11) 
with the control chosen according to the third of equations (1.12) 
will produce the optimum path. At some t=T on this path, the con- 
straints will be satisfied, the desired variable will be maximiz- 
ed, and simultaneously (1.12) and (1.14) will be satisfied. 

In Chapter III a method will be given that will allow us to 
correct from an arbitrary point U 1 ( 0 ) toward U* ( 0 ) provided U ' ( 0 ) 
is within some unknown region R in the vicinity of U*(0). The 
difficulty is to find some point U*(0) which lies within R. The so- 
called Direct Method of solution used here reduces the problem to 
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the solution of a two-point boundary-value problem involving the 
extremals. Consequently, the missing boundary conditions must 
somehow be guessed. In general, the guess must be within some 
region R, if the corrective program is to converge. In the Grad- 
ient, or Steepest Ascent methods of Kelley and Bryson a nominal 
trajectory is guessed and the two-point boundary-value problem 
is solved by correcting along paths which are not extremals. This 
method too has its disadvantages. For some problems the optimum 
path may never be reached, although the correction program "con- 
verges", that is, satisfies the convergence criterion. For the 
problem treated in this thesis, a path was obtained by the direct 
methods given here which attained a final value of the variables 
to be maximized well beyond that of the "optimum" path obtained 
by the method of Steepest Ascent. 

2 . 2 Nominal paths . 

An advantage of the Gradient Method is the following. In 
physical optimization problems vie often have some idea beforehand 
what kind of control program is likely to come fairly close to the 
optimum one. The only "guessing" involved in the Gradient Method 
is to choose such a program and construct a nominal path which 
connects the given initial constraints with the given terminal 
constraints. Provided the guess is good enougn, the path is cor- 
rected by making the integrated value of H smaller and smaller 
along the path. 
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In a given problem, the degree of "goodness" required of this 
nominal path when using the Gradient Method corresponds in a rough 
way to the "goodness" of the first guess for U'(0) when using extre~ 
mals. In the former case, however, one can be guided by physical 
reasoning. The physical meaning of the value of U* ( 0 ) is by no 
means as clear. 

In the preliminary work associated with this thesis, the authors 
have worked on a means to generate a point U * ( 0 ) within R from a 
nominal path. In the sample problem undertaken, the method works 
well and provides a value for U'(0) which is quite close to U*. 

The method is quite simple and straightforward and adds very little 
complexity to the computer program which is used to correct from 
U»(0) to U*(0). 

2.3 The fundamental adjoint set 

As a primary tool, we make use of the fundamental set which 
is used by Faulkner ( 5 )• The fundamental set is the nxn matrix, 



( 2 . 1 ) 



U = 



U 11 u 12 



U 21 U 22 



• • 



u nl u n2 



u 



In 



0 . u 



2 n 



. . u 



nn 



whose value at t=0 is the nxn identity matrix and whose elements 
are chosen so that u. . is the j'th component of the i'th solution 
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to the adjoint equation (1.5). Note that each row vector of the 
fundamental set is a particular solution to the adjoint (1.5). 
For instance, 



( 2 . 2 ) 




Now, consider a particular combination of the row vectors U^, 

(2.3) U = c x U 1 + c 2 U 2 + . . . + c n U n , 

where the c's are constants. 

This may be written as 



U = CU, where 

C = [Jl C 2 * * * C n}‘ 

Note that as thus defined, U(0) = C, so that the C's are the initial 
values of the adjoint variables. At any time t, u^ is the ith com- 
ponent of U: 



(2.4) 

(2.5) 



( 2 . 6 ) 

or 



u i = °1 u li + °2 u 2i + 




u 



nx 



(2.7) \ = Cjx , Where 

u. is the ith column of the fundamental set. 

“X 

One of the c's may be chosen arbitrarily. If some state variable 
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does not appear explicitly in the equations of motion, then the 
corresponding adjoint variable is constant. If this state vari- 
able is to be maximized then it will be convenient to make the 
corresponding c equal to plus one. If the variable is free at 
either end point then the corresponding adjoint variable is zero 
everywhere, and a different c will be chosen arbitrarily. In any 
case, we must be sure that the corresponding adjoint variable is 
non-zero . 

Program for generating the C- vector from a nominal path 

In this section a program is given for generating a starting 
set of constants U'. A nominal control program is chosen from 
physical reasoning which appears likely to be fairly close to 
optimal and which will generate a path from the given initial 
conditions to a terminal point which satisfies or nearly satis- 
fies the terminal constraints. This may be simple or difficult, 
depending on the problem. If such a program cannot be found, 
then this method for starting cannot be used, nor can be Gradient 
Method, since both are identical up to this point. In this case, 
we must go back to the "guessing game". It may be necessary to 
"map" the U(0) hyperplane to find which regions give good start- 
ing points, and then to try various points within those regions. 
That is, guess the c*s and calculate the corresponding trajec- 
tories in some systematic way. 



25 



Suppose however that a satisfactory nominal control program 
is found.. The equations of motion are now integrated numerically 
from the given initial conditions using the nominal control pro- 
gram, stopping when a terminal constraint is met. At the same 

p 

time, we integrate also the n equations of the fundamental set. 

At the terminal point S(T), each of the elements of the fundamental 
set has some value u. .(T). 

If the path had been the optimum path, equations (1.14) 
would have been satisfied at the point S(T). We will choose the 
vector C so that they are satisfied, using the values from the nom- 
inal path. Then this vector C becomes U'(0) for starting the cor- 
rective iterations. If the program does not converge, the nominal 
path may be varied somewhat to get another trial C„ In this way, 
the problem of guessing the initial values of the adjoint variables 
may be eliminated, as in the Gradient Method and still the optimum 
path is approached through extremals. 

It is to be noted that equations (1.14) are not the only 
relationships which are used to furnish the required number of 
equations in the C's. According to the particular problem, the 
Hamiltonian may be constant, or a particular adjoint variable may 
be constant. These allow use of relations calculated for t=0 
rather than at t=T as above. 
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2.5 Faulkner's method of differential corrections. 

A method was given above for obtaining a first approximation 
for the starting values of the adjoint variables from a nominal 
path. Let us now suppose that such a set is in hand which was 
obtained in this way or in some other way. 

If the equations of motion (l.l) together with the adjoint 
equations (1.5) are now integrated numerically with the control 
chosen according to the maximum principle or (1.7) a path results 
which is an extremal in that it furnishes an optimum to some end 
state. But the path obtained using this particular set of initial 
conditions does not, in general, satisfy the terminal constraints 
(1.2) and the transversal conditions (1.14). We must have some 
means for correcting these constants (which we have variously called 
the C-vector or the vector U(0) in earlier sections) so that (1.2) 
and (1.14) are satisfied at some point S(T)„ 

The method used here is the method of differential corrections 
developed by Faulkner (ref. 5) which makes use of differentials in 
an iterative scheme which is of the N ewton-Raphs on type. Making use 
of the fundamental set and the development of Chapter I, the method 
may be presented quite simply. 

Equation (1.6) provides the relation (for a problem without cor- 
ners) 




UF <$Pdt 
P 
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For every state variable s i which is specified at t = 0, &s^(0) = 0. 
For those which are free, we will require u^(0) = 0, so that u^<5sj. = 0 
for all cases. This procedure takes care of separated end conditions 
in a simple and straightforward manner. This case causes a great deal 
of difficulty with some methods of solution. Note that there are 
still n unknowns at the start, since the value of s^ is unknown but 
the value of the corresponding u^ is known. 

With this simplification, we have 

(2.8) 6 <5s] t = fj (UF p ) 6P dt, 

where the dimensions are nxn for U, nxl for &S, nxm for F p and mxl 
for P. From (1.7) we will choose the p’s so that 
H p = UFp = 0. 

Now if S and t are fixed and the c's are varied by small amounts &C, 
then P must vary so that 

(2.9) 6cU c F p + 6P T H pp = 0, 

where the superscript T indicates the transposed matrix. But since 
U = CU , then U c = U. Making this substitution, and noting also 
that H is the symmetrical mxm matrix whose ij'th element is 

irir 

< Vij = 3%/5p i 3p J ’ 

we now form the mxl transposed of (2.9) • 

(2.10) (UF p ) T <5c T + H pp 6P = 0 
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This equation is solved for &P, and the result is substituted 



into (2.8). 



6p = - (H ) _1 (UF ) T <$c T 

pp — p 



( 2 . 11 ) 



T —ITT 

[u Ss} = - J (UF p ) (Hpp) (UF p ) Sc dt, 
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So 
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_ _J 



,T -1 T . T 

It, 

T 'o 

a form convenient for machine programming. 

Novi suppose each term in the integral on the right be denoted 

as I .. . lie have then at t = T 
J 



( 2 . 12 ). 



Since one c is arbitrary, say c-^, then <5c^ = 0 and the equation 
furnishes n-1 relations among the 6c' s and the <£s's. 

In addition, the transversal conditions provide n-k equations 
whose variational form is, for the problem to be considered here, 

£u^ = 60^. 

Now for corrections we let 

(2.13) ~ - u i(T) - ^(T) 6 t = 

and also, for the k constrained variables, 

6s ± = - Si (T) - s^T) 6 t + s c [T) 
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where s (T) is the given terminal constraint for the ith variable. 

c i 

Since there are k of these variables constrained, there remain 
(n- k_l) elements &s^ in equation (2.12). These may be eliminated, 
yielding k equations. The n-k equations (2.13) are adjoined to 
these, for a total of n equations involving the unknowns 6 t, 6 c ?, 

6c^, . . . 5c n in terms of the integral elements I^» the 
elements of U(T), the deriviatives at t=T and the differences (2.13). 
The equations are solved for the unknowns, and these are applied as 
corrections for the start of the next iteration. 
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Chapter III 

DEVELOPMENT OF THE PROBLEM 

3.1 General description of the problem and assumptions. 

The problem described in this thesis is that of deter- 
mining the optimal angle of attack program p(t), which will 
maximize the distance covered over the earth's surface by a 
hypersonic glider or lifting vehicle which has been injected 
into some initial re-entry point. The starting point of the 
trajectory is specified by the re-entry injection parameters, 
the initial conditions, which are: 

v(0) = 35,920 feet/second - initial velocity 

h(0) = 30^,000 feet - initial altitude 

(3-D 

x(0) = 0 nautical miles - initial distance 

z(0) = 0.18 degrees - initial flight path angle 

These parameters were taken from reference (8) so that a 
comparison could be made between the Gradient Method used by 
earlier authors and the method of differential corrections 
which was used in this thesis. Also, the value of wingload- 
ing, the value of acceleration due to gravity, the value for 
a standard earth radius, and the value of a standard nautical 
mile were taken from the same reference. 

V’ingloading = mg n = 27. 3 lb/ft^ , where: 

A 

m = mass of the vehicle 
A = wing plan form area 

gQ = 32.2 ft/sec - acceleration due to gravity at 
the earth's surface 

The model atmosphere used was based on the tables of reference 
(9). ^he atmospheric density was assumed to have the form 
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• . y 

, v.'here the parameter b was calculated to fit the tab- 
ular values at h = 0 feet and h = 200,000 feet. It was fur- 
ther assumed that the atmosphere was spherically symmetric and 
fixed with respect to the earth for simplifying reasons. 

3.2 Equations of motion. 

The coordinate system used is depicted on Fig.3.1and the 
equations of motion are: 

s = Eh ’ cos 2 ’ 



(3.2) 



ft = v sin z , 



D 



V = - - - g cos z , 
m ° 



L v 

z " tnv + R+h 



cos z 



_ £ 



cos z 



w nere : 



R = 3440 nautical miles, the radius of the earthy 
d = C L ) , 

L = |(/»v ? A C D ) , 
r R -;2 

S = SoLRThJ * 



The coefficients of lift and drag are from reference (8) and 
are : 

C D = C LL l sin 'pl + C do 

Where: C L0 = 1.82 ; C DL = 1.46 ; C DC = 0.042 ; and p is the 

angle of attack of the vehicle. 

3.3 Adjoint equations. 

r Lhe system of differential equations which is adjoint 
to the variational equations of the equations of motion (3.2) is 
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The fundamental set of adjoint equations U or equivalently 
(u^) was chosen such that the fundamental set at time t = 0 
is the identity matrix of rank four. That is: 



(3.4) 



2 = (u ij> = 
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U zJ 



( ■n r p o 

^ q =7q (u,.)dt , where the u^ are defined 
by the adjoint equations (3*3) * We have the further relation 
that: U = [u y u^ u y u z J = [c-, c ^ c^ c^][u]*^ which is a more 
convenient form for computations using Fortran programming, 

Tn general, the solution of the adjoint equations is 
only determined to within a multiplicative constant , hence 
we can choose one of the c's. Tn this problem, we chose c^=l. 
If we consider the equations of motion in vector or 
matrix form, we have the following relationships: 
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s = 



x 

h 

v 

_z J 



P = S ; F p= 



0 

0 

“Dp 

m 

Dp 

mv, 



-i. f r> 

where = ~=r- 
P cp 



The Hamiltonian is given by H = U*F. 

J.b Maximum principle. 

For any extremal, we must satisfy the condition that 
U.Pp = o for any point on the extremal. This condition 
uniquely determines the angle of attack program for that 
particular extremal since: 

IK "4§ u v - o 

Substituting for 4^ and 42 and solving for p we obtain: 



(3.5) 



p = arc tan 
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C u z 
°DL Z 



T 



- u. 



f °r - f < P < 2 



8.5 Pnd conditions. 

Since we wish to maximize the total distance traveled 
over the earth's surface, we must have the following conditions 
at the end time t = T: 
x(T) = maximum 
(3.6) h(T) = 0 

v(T) , z(T), and T free. 

3.6 Variational equations. 

On any curve without corners, we have the following 

relation between the adjoint and the variations: 

T 

(3»7) [u^6x + u^6h + u^6v + u^6z] t =^ fcp dt j =1,2,3,^ 



For an extremal we must also satisfy the Euler equations or 
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the maximum principle of Pontryagin so that we have 

r r i'-» 

l‘U*P dt = maximum and 0 along the extremal. If we 

*T) 

consider x,h,v,z, and t fixed and change the constants then 



PU 
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+ TT 6C, 
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SCr> 2 0c o ww 3 ^ 6c, , h 1 ~p ' " " pp 



.] 



M • F + U*P 6p = 0. Noting that 



= U" , i =2,3,4, we can solve for 6p and substitute in 
equation (3.7) which will give us four equations with 

which we can correct the initial values of the c's for a 

terminal variation of the altitude h: 

[ 1 ^ T j i 

u^6x + u^6h + u^6v + u”,6zj -3~ J (U*F p )(U*F p ) 

i=2 r\ z ^ ~r “ 
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Since u = 1 , u = u' = u = 0 for all t, and noting that 

X XXX 

*ij = we can solve the following set of equations for 6h 

in terms of the 6c, : 
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Prom the transversal conditions, we have that: 
u y (T) = 0 
(3.10) u (T) = 0 

(xu ? +hu h ) T = 0 

Prom the last equation of (3.10) we have: u^(T) = (-cot z). 



Since the Hamiltonian H = U*F and H(T) = 0 from equations 
(3.10), then the Hamiltonian is identically zero for all t. 
This fact will be used to determine one of the constants (c 2 ) 
and effectively reduce the nroblem of guessing the initial 



c's to start the calculations. 
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The transversal conditions (3*10) give us three more 
independent variational equations associated with the end 
conditions : 

fiu v (T) = [u^£c 2 + u^6c 3 + Uy£c i .] T 
(3*11) fu z (T) = [u z £c 2 + u^6c 3 + u z 6c 4 l T 

6[xu x + hu h J T = h(T)[u^6c 2 + u^5c 3 4 u h f ' c ijJ T 

3.7 Differential corrections 

The method of differential corrections given by 
Faulkner in reference (5) was used for the solution of the 
optimization oroblen. In general, we will guess the initial 
values associated with the adjoint enuations and the terminal 
time T. For this problem, we must choose eg * c 3 , c 4 , nd T. 

We then correct these parameters by calculating an extremal 
with the initial guess and then solve for corrections to the 
parameters in terms of the end values associated with the 
extremal. For all end conditions of the form s(T) = S, we will 
generally have some error in s(T), We set 6s (T) = S - s('T) - 
s(T) 6T. For this problem, we set 6h(T) = h(T)-h(T)6T and in 
a similar manner, we obtain from equation (3*11 ) ° 

-u v (T) = [ujec 2 + ujtc 3 + u^6c 4 ] t + u v (T)6T 

(3*12) -u z (T) = [u^c 2 + u^£c 3 + u z 6c 4 ] t + u z (T)6T 

* 2 3 j[i 

_COt Z J ,p = h ( T ) w U^t u h^ °3 d" ^ J rp 

+Cxu 4 hu h 4 hu h 3 T 6T 

These equations together with the solution to (3*9) for 6h in 

terms of the 6c' s will give a set of equations from which we 
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solve for the fc's and fT to make corrections in the starting 
values for the subsequent trajectory. 

3.8 Computation of approximate starting values by use of 



nominal trajectories. 

In general, for a nth order set of differential 
equations, one must choose (n-1) of the constants of the sol- 
ution to the adjoint equations to start the problem. This is 
generally a difficult problem in itself, since the constants 
are not related to the ohysical aspects of the problem in any 
discernable manner. One of the primary aspects of this thesis 
was to find a method which could be used to find a logical 
choice for the initial valves of these constants. Fortunately 
the Hamiltonian was identically zero in this problem. This 
fact, together with the transversal conditions (3«10) enabled 
us to solve for approximate constants by first computing a 
nominal trajectory; here we choose a nominal angle of attack 
program p(t), which seems likely in a physical sense to 
give a maximum distance trajectory. The first such nominal angl 
of attack program tried was to use the maximum lift over drag 
ratio where p=20.5 degrees, a constant. This did not work. 

We then used the following angle of attack program for which 
we are indebted to Professor Faulkners 



p(t) = 5^.7 



1 - h(t) 



750,000 



The value of 5^*7° for P is the value at which we have 
maxium lift. It was supposed that the angle of attack was 
never negative and from the data of constant angle of attack 
programs, the value of 75O,OO0feet was chosen to insure that 
p was never negative. 
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solve for the fc's and t T to make connections In the starting 
values for the subsequent trajectory, 

3.8 Computation of approximate starting values by use of 
nominal trajectories. 

In general , for a nth order set of differential 
equations, one must choose (n-1) of the constants of the sol- 
ution to the adjoint equations to start the oroblera. This is 
generally a difficult problem in itself, since the constants 
are not related to the ohysical aspects of the problem in any 
discernable manner. One of the primary aspects of this thesis 
was to find a method which could be used to find a logical 
choice for the initial valves of these constants. Fortunately 
the Hamiltonian was identically zero in this problem. This 
fact, together with the transversal conditions (3«1°) enabled 
us to solve for approximate constants by first computing a 
nominal trajectory; here we choose a nominal angle of attack 
program p(t), which seems likely in a physical sense to 
give a maximum distance trajectory. The first such nominal angle 
of attack program tried was to use the maximum lift over drag 
ratio where p=20.5 degrees, a constant. This did not work. 

We then used the following angle of attack program for which 
we are indebted to Professor Faulkner" 

p(t) = 5^.7° [l - h(t) 1 

l 750,000 J 

The value of 5^*7° for P is the value at which we hove 
maxium lift. It was supposed that the angle of attack was 
never negative and from the data of constant angle of attack 
programs, the value of 75°»°°0feet was chosen to insure that 
p was never negative. 
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With this choice for the angle of attack program, we 
obtained a good set of starting constants by integrating the 
fundamental set U along with the equations of motion until h = 0, 
At this time, we solve the set of equations: 
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for C 2 , , and 04 . These equations are the transversal 

equations (3 . 10) . This problem of guessing the c 8 s is further 

reduced due to the Hamiltonian being identically zero, since 

for a given choice of c^ and 04 , the angle of attack is uniquely 

determined from equation ( 3 * 5 ) and then we may solve for 02 from 

the Hamiltonian at time t = 0. If we have the case where h(0)= 

0 , then we have a singularity and connot solve for C£ , but in 

this case, the choice of 02 is arbitrary. Effectively, this 

problem is then reduced by an order of one and we need only 

guess two constants. We felt that if we had an extremal, we 

would determine its constants in this manner and hence we could 

which 

get estimates for the constants from a tra jectoryj^was a "good" 
approximation to an extremal. 

3.9 Computational scheme. 

This problem was orogrammed on a CDC 1604 computer 
using a Runge-Kutta integration method. The computational 
scheme is as follows: 

i) We need a nominal trajectory to start. To get it, we 
first intergrate the equations of motion 0 . 2 ) and the fundamental 
set ( 3 . 4 ) using a programmed angle of attack and ending the 
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trajectory when h = 0. At this time, we solve for the 
initial constants 03 , , and c^ from equation ( 3 . 13 )* We then 

”se these constants to start the first optimal trajectory. 

ii) Then we compute the first optimal trajectory where 

the angle of attack program is determined by the maximum 
principle. We integrate the equations of motion, the fundament** 
al set, and the nine in equation ( 3 . 9 ) a total of twenty- 

five equations. 'We terminated the integration when one of 

the following conditions is met: h = 0 , u v = 0 , u^ = 0 , or 

I z I = g, since it was felt that the end of the significant part 
of the trajectory had been reached. We felt that in this 
problem no optimum would occur for u^< 0 since the resulting 
high drag would dissioate energy needlessly. The supposition 
that this is a maximum distance trajectory dictates that lzl< 

Since we do not have t appearing explicitly in the equat- 
ions and since the terminal conditions determine T, we do not 
use 6 T in our calculations. We only need the 6 c ' s the correct 
the c's to start the next trajectory. 

iii) Vie then repeat the computations as in the above 
paragraph, using the new c's for the next trajectory. Vie 
test for convergence at the end of each trajectory: 

(3.1*0 If [h(T)j 2 + [xu y + hu h ] 2 + [u v (T )] 2 + [u z (T )] 2 =£, 

we say that we have converged to the solution. 

3.10 Results and conclusions. 

The optimum trajectory that was calculated is deDicted 
on Pig. 3.3 and the corresponding angle of attack program is 
denicted on Fig. 3 A On these two figures, a comparison has 
been made between the results of this thesis and results from 
reference ( 8 ) , which were obtained by the Gradient (or Steeo- 
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est Ascent) Method. The comparison can r.ot be considered as 
an exact one since the integration method, the time step 
ofc integration interval, and the approximation to the 
atmospheric models may not be identical. However, there 
appears to be a significant increase in the maximum distance 
obtained using the method of differential corrections. 

The convergence criteron ( 3 . 14) we initially chose turned 
out to be a poor choice. We could not determine how close we 
were to the optimal solution. We then proceeded to "map' 1 
the plane of C 3 and 04 by using a systematic choice of 03 and 

c^ and computing the corresponding extremals. The results of 
the mapping procedure are given on Pig. 3.5, Fig. 3. , and Fig„3c7° 
Fig. 3.5 shows the contours of distances obtained by extremals 
as a function of c.^ and c^. Fig. 3*6 gives detailed contours 
around the optimal solution. Fig. 3.7 is a cross section of the 
"ridge” which occurs in the Cj and c^ plane. It is a con- 
jecture of the authors that some difficulties are encountered 
by a gradient method solution when there exists such a ridge 
in a problem, since the slope is nearly zero. 



The major problem encountered was the apparent instability 
of the problem near the end of an extremal, W e attempted to 
integrate backwards from the terminal conditons and found that 
this was virtually impossible. The problem seems to be unstable 
we could not even integrate backwards with constant p. Back- 
ward integration was very good if we started backward from a 
point that was not too close to the end of the trajectory. 

Another main problem encountered was that of determining 



40 



the end time on a trajectory. We first thought that all we 
need consider is the time when h = 0 . However, we found 
tha- other conditions effectively terminated the useful part 
of the extremal earlier. For example, most extremals led to 
a point where |z|> —• 

Theoretically, on the maximum distance trajectory, all 
terminal conditions would be met simultaneously. In the actual 
computations, this did not occur. The most probable reason 
that we have this apparent discrepancy is that the round off 
errors in the calculations and the inherent error in the 
integration method are large enough to prevent us from satis- 
fying all conditions at the end time. 

In the final analysis, we found that a program which 
reduced the magnitude of the corrections would lead to con- 
vergence. The method used was to calculate a set of c's 
from a nominal trajectory and then calculate the first tra- 
jectory and the fc's corresponding to that trajectory. V e 
store x, the c's, and 6c* s of this trajectory and then cal- 
culate the next trajectory. If the distance of this trajectory 
is greater than the distance of the preceding one, we proceed 
with the iterative process. If the distance is less, then 
we reduce the corrections by one half and calculate the tra- 
jectory again. If we still have less distance, we further 
reduce the corrections until we obtain a greater distance. 

The reason we developed this method was because the magnitudes 
of the fc’s were of the order of magnitudes of the c's 
was noted that the integrals relating 6h to the fc's are 
improper integrals. This, in addition to the unstable char- 
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acter of the end of the trajectory, tended to make the 
magnitude of the 6c' s large; however, the program did 
correct in the right direction. These integrals seem to be 
divergent for the extremals which satisfy the transversal con 
ditions. A formulation of the oroblem suggested by Professor 
Faulkner which would avoid these divergent integrals by using 
a set of u's which were determined at the end time was pro- 
grammed. However, this program was not checked out or run 
since the original integrals did not actually diverge in cal- 
culations made on the computer. The method suggested was to 

■K- * 

use a fundamental set U such that U (T) = . After com- 

1 J O 

puting an extremal using U , we solve for the matrix of con- 
stants [a*^] which satifies the relation: U(T) = ^a^ ^ ]U(T) . 

- A * 

We then compute the same extremal using U = U. + cU„ , where 

Ll l £ 

-A \ 

U. = 2_ a - -U , i= 1,2. Ve then have the relation 

. 1=1 1J 



£h(T) 




m k 

( u 0 • F ) 
' 2 o 



U • F 



PP 



2 

dt be 



and we solve for be to correct the 



starting vaiues. Further, this integral is not divergent. 

The first integration method used in calculations was 
rectangular integration and it was found that trajectories 
thus calculated were very rough. A comparison of the rectang 
ular integration method and the Runge-Kutta method was made 
and we chose the latter method due to its accuracy and 



smoothness . 
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Figure 3.1 




Coordinate System 
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Figure 3*3 * ^- T 
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Figure 3*^ m. 
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Optimum Angle of Attack Program 



Figure 3-5 
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ontours are in thousands of nautical 
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Figure J.6 



Detailed Contours of the c 0 - plane 



:ontours are given in thousands of nautical mil: 
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. Figure 3» 7 
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